# This is a sample Python script.

# Press Shift+F10 to execute it or replace it with your code.
# Press Double Shift to search everywhere for classes, files, tool windows, actions, and settings.
# See PyCharm help at https://www.jetbrains.com/help/pycharm/
#coding=utf-8
import random
from urllib.request import urlopen
from urllib.request import Request
import numpy as np
import json
import geopandas as gpd
from shapely import geometry
import matplotlib.pyplot as plt
from shapely.geometry import Polygon as geoPolygon;
from geopandas import GeoSeries, GeoDataFrame
from math import sin, asin, cos, radians, fabs, sqrt
EARTH_RADIUS = 6371

##GeoJSON = json.loads('{"type":"MultiPolygon", "coordinates":[]}')
GeoJSON = json.loads('{"type":"FeatureCollection", "features":[]}')
#LineJSON = json.loads('{"type":"Feature", "geometry":{"type": "LineString","coordinates": []}}')

#Linecoordinates = [] #保存线
Polycoordinates = [] #保存面
boundsAreaList = [] #保存边界面积
numStat90 = 0
numStat80 = 0
numStat50 = 0
numStat20 = 0

def hav(theta):
    s = sin(theta / 2)
    return s * s

def save_json(save_path,data):
    assert save_path.split('.')[-1] == 'json'
    with open(save_path,'w') as file:
        json.dump(data,file)

def get_distance_hav(lat0, lng0, lat1, lng1):
    "用haversine公式计算球面两点间的距离。"
    # 经纬度转换成弧度
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lng0 = radians(lng0)
    lng1 = radians(lng1)
    dlng = fabs(lng0 - lng1)
    dlat = fabs(lat0 - lat1)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))

    return distance

def readBoundsJSON():
    index = 1
    for bound in open("land_bound.txt"):
        jsonObj = json.loads(bound)
        boundsObj = jsonObj["bounds"]
        if ( len(boundsObj) > 0 ) :
            boundsArea = jsonObj["bounds_area_size"]
            userId = jsonObj["user_uid"]
            pointsList = boundsObj[0]["points"]
            polygon = []  # 保存一个边界数组
            member = []  # 保存一个点数组
            Polycoordinates = []
            for point in pointsList:
                pointlanLat = str(point["lng"]) + "," + str(point["lat"])
                member.append(pointlanLat)
                polygon.append(member)
                member = []

            print(index)
            index = index + 1

        tempPolyJSON = json.loads('{"type":"Feature", "geometry":{"type": "Polygon","coordinates": []},'
                      '"properties": {"boundsArea":"","userId": ""} }')
        if( index == 2 ):
            tempPolyJSON["properties"]["boundsArea"] = boundsArea
            tempPolyJSON["properties"]["userId"] = userId
            tempPolyJSON["geometry"]["coordinates"].append(polygon)
            GeoJSON["features"].append(tempPolyJSON)
            boundsAreaList.append(boundsArea)
        else:
            tempGeoJSON = json.loads('{"type":"FeatureCollection", "features":[]}')
            tempPolyJSON["properties"]["boundsArea"] = boundsArea
            tempPolyJSON["properties"]["userId"] = userId
            tempPolyJSON["geometry"]["coordinates"].append(polygon)
            tempGeoJSON["features"].append(tempPolyJSON)
            if(True): #SpatialAnalysis(GeoJSON, tempGeoJSON)
                GeoJSON["features"].append(tempPolyJSON)
                boundsAreaList.append(boundsArea)

def replaceGeoJSONStr(GeoJSON):
    GeoJSONStr = json.dumps(GeoJSON)
    GeoJSONStr.replace('["','[')
    GeoJSONStr.replace('"]', ']')
    GeoJSONResult = json.loads(GeoJSONStr)
    return GeoJSONResult

def SpatialAnalysisTest():
    p1 = geoPolygon([(0, 0), (1, 0), (1, 1), (0, 1)])
    p3 = geoPolygon([(0.9, 0), (3, 0), (3, 1), (1.0, 1)])
    p2 = geoPolygon([(2, 0), (3, 0), (3, 1), (4, 1)])
    array = []
    array.append(p1)
    array.append(p3)
    geo1 = gpd.GeoSeries(array)
    geo2 = gpd.GeoSeries(p2)
    isCrosses = geo1.disjoint(geo2)  # 判断相交

    fig, ax = plt.subplots()
    geo1.plot(ax=ax, color='orange', markersize=100, marker='*')
    geo2.plot(ax=ax, color='orange', markersize=100, marker='*')
    plt.xticks(rotation=20)
    plt.show()

    print(isCrosses)

def isEqualsTwoPolygon(GeoJSON, tempGeoJSON):
    features = GeoJSON["features"]
    tempfeatures = tempGeoJSON["features"]
    tempAllPolygons = []
    nextPolygon = []
    nextPolyArea = 0
    isAddPolygon = True
    for tempfeature in tempfeatures:
        poly = tempfeature["geometry"]["coordinates"][0]
        nextPolyArea = tempfeature["properties"]["boundsArea"]

    # 生成所有多边形的geoPolygon
    for feature in features:
        poly = feature["geometry"]["coordinates"][0]
        fristPolyArea = feature["properties"]["boundsArea"]
        if (fristPolyArea == nextPolyArea):
            isAddPolygon == False
            print("isEquals")
            return isAddPolygon

    return isAddPolygon

def SpatialAnalysis(GeoJSON, tempGeoJSON):
    features = GeoJSON["features"]
    tempfeatures = tempGeoJSON["features"]
    tempAllPolygons = []
    nextPolygon= []
    nextPolyArea = 0
    nextpointX = 0
    nextpointY = 0
    firstpointX = 0
    firstpointY = 0
    isAddPolygon = True
    for tempfeature in tempfeatures:
        poly = tempfeature["geometry"]["coordinates"][0]
        nextPolyArea = tempfeature["properties"]["boundsArea"]
        polyStr = []
        for point in poly:
            nextpointX = float(point[0].split(',')[0])
            nextpointY = float(point[0].split(',')[1])
            polyStr.append(point[0].split(',')[0])
            polyStr.append(point[0].split(',')[1])

        reshapeNum = int(len(polyStr) / 2)
        nextPolygon = geoPolygon(np.array(polyStr).reshape(reshapeNum, 2))

    geoNext = gpd.GeoSeries(nextPolygon)

    # 生成所有多边形的geoPolygon
    for feature in features:
        poly = feature["geometry"]["coordinates"][0]
        fristPolyArea = feature["properties"]["boundsArea"]

        polyStr = []
        for point in poly:
            firstpointX = float(point[0].split(',')[0])
            firstpointY = float(point[0].split(',')[1])
            polyStr.append(point[0].split(',')[0])
            polyStr.append(point[0].split(',')[1])

        #reshapeNum = int(len(polyStr)/2)

        tempPolygon = geoPolygon(np.array(polyStr).reshape(int(len(polyStr)/2), 2))
        geofirst = gpd.GeoSeries(tempPolygon)
        isAddPolygon = TwoPolygonOverlay(geofirst, geoNext)
        if(isAddPolygon == False):
            return isAddPolygon
        # tempAllPolygons.append(tempPolygon)
    # geoAll = gpd.GeoSeries(tempAllPolygons)
    # fig, ax = plt.subplots()
    # geoAll.plot(ax=ax, color='orange', markersize=100, marker='*')
    # geoNext.plot(ax=ax, color='red', markersize=100, marker='*')
    # plt.xticks(rotation=20)
    # plt.show()

    return isAddPolygon

def CompareTwoPolygon(geofirst, geoNext):
    # geofirstConter = geofirst.centroid
    # geoNextConter = geoNext.centroid
    # pointDist = get_distance_hav(float(geofirstConter.x), float(geofirstConter.y), float(geoNextConter.x),
    #                             float(geoNextConter.y))
    isAddPolygon = True
    # isintersects = geofirst.crosses(geoNext)  # 判断相交
    isEqualses = geofirst.equals(geoNext)  # 判断相等
    # isintersects = geofirst.intersects(geoNext)  # 判断相交
    if (isEqualses == True):
           isAddPolygon = False
           print("isEquals")

    return isAddPolygon

def TwoPolygonOverlay(geofirst, geoNext):
    minArea = min([float(geofirst.area), float(geoNext.area)])
    df1 = gpd.GeoDataFrame({'geometry': geofirst, 'df1': [1]})
    df2 = gpd.GeoDataFrame({'geometry': geoNext, 'df2': [1]})
    res_intersection = gpd.overlay(df1, df2, how='intersection')
    isAddPolygon = True
    if(res_intersection.empty != True):
        if( float(res_intersection.area) > float(minArea)*0.9):
            print("相交大于90%")
            isAddPolygon = False
        elif (float(res_intersection.area) > float(minArea)*0.7):
            print("相交大于70%")
            #isAddPolygon = False
        elif (float(res_intersection.area) > float(minArea)*0.3):
            print("相交大于30%")

    return isAddPolygon


if __name__ == '__main__':
    ##print_hi('PyCharm')
    readBoundsJSON();
    save_json("./GeoJSON.json", GeoJSON)
    # GeoJSON = replaceGeoJSONStr(GeoJSON)
    np.savetxt("./boundsAreaList.txt", boundsAreaList, fmt='%.08f')
    # GeoJSONData = gpd.read_file("./GeoJSON.json")
    print("over")
    # SpatialAnalysis(GeoJSONData)
    # SpatialAnalysisTest()

